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ABSTRACT 


A  technique  for  reducing  the  number  of  iterations  necessary 
for  solving  linear  programs  using  the  primal-dual  algorithm  is 
presented.   It  appears  that  the  new  method  will  also  decrease  the 
number  of  iterations  over  any  other  simplex  algorithm.   A  FORTRAN 

program  incorporating  the  technique,  as  well  as  some  comparative 

v 
computational  results  are  given. 
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SECTION  I 
INTRODUCTION 

The  primal-dual  algorithm,  as  developed  by  Ford  and  Fulkerson, 
theoretically  permitted  faster  convergence  in  solving  linear  programs 
over  such  methods  as  the  two-phase  or  Charnes '  M-method.   Instead  of 
first  driving  the  artificial  variables  to  zero  to  obtain  a  feasible 
solution  and  then  working  to  optimality  in  a  second  phase,  the  primal- 
dual  works  towards  optimality  and  feasibility  in  the  primal  and  dual 
problems  at  the  same  time.   This  is  accomplished  during  what  would 
be  the  first  phase  for  the  other  procedures.   However,  in  actual 
computational  experience,  the  primal-dual  algorithm  has  not  achieved 
faster  convergence  over  the  other  simplex  procedures. 

This  paper  presents   a  modification  of  the  primal-dual  algorithm 
that  does  accomplish  faster  convergence  than  the  usual  primal-dual 
algorithm  and  perhaps  also  the  other  simplex  methods. 


Professor  Harold  Greenberg  provided  the  idea  for  the  new 


technique. 


SECTION  II 


THE  STANDARD  PRIMAL- DUAL  ALGORITHM 


Using  matrix  notation,  the  linear  programming  problem  is: 

Minimize  z  =  z  +  CX 
o 

Given    AX  =  H 

x.^0     j=l,. . . ,n 
J 

where  z   is  a  constant  and 
o 

C  is  a  1  x  n  vector  with  components  c.  ^  0,  j  =  l,...,n 

H  is  an  m  x  1  vector  with  components  h.  2:  0,  i  =  l,...,m 

A  is  an  m  x  n  matrix  with  elements  a.  .   i  =  l,...,m 

i  >  J  > 

j  =  1,..  .,n 

X  is  an  n  x  1  vector  with  components  x.,  j  =  l,...,n 

2 

The  standard  primal-dual,  as  it  appears  in  Dantzig,   consists 

in  augmenting  the  constraint  equations  with  artificial  variables  as, 


AX  +  Y  =  H  (1) 

y.^0         i=l,...,m 


where  Y  is  an  m  x  1  vector  with  components  consisting  of  the  arti- 
ficial variables  y.;  and  forming 


w  =  w  +  DX 
o 

where 

m 
w  =  S  h . 
°   1-1  X 

and  D  is  a  1  x  n  vector  with  components  d.,  j  ■  l,...,n  where 


2 
George  B.  Dantzig,  Linear  Programming  and  Extensions 

(Princeton:   Princeton  University  Press,  1963),  pp.  247-53. 


m 
d.  =  -  E  a. 
J    i=1   i.J 

It  is  noted  that  both  the  c.  and  the  h.  must  be  positive 

J  * 

prior  to  starting  the  computations.  Negative  h.  may  be  made 
positive  by  multiplying  the  associated  constraint  equations  by 
a  minus  one.   The  method  for  correcting  negative  c.  is  discussed 
in  the  "Description  of  the  Computer  Program"  section  of  this 
paper. 

The  problem  is  to  find  x.  ^  0,  w  =  0,  that  minimizes  z 
subject  to  (1).   At  any  stage  of  the  iterative  procedure  the 
equations  have  the  form: 

XB  +  AXX  =  Hx 

Y  +  A£X  =  H2  (2) 

DX  =  w  -  w 
o 

CX  =  z  -  z 
o 

where  X   is  a  vector  representing  the  current  basic  variables 

D 

and  X  is  a  vector  representing  the  current  non-basic  variables. 
The  numerical  subscripts  1  and  2  represent  partitioned  matrices 
and  vectors.   Matrices  and  vectors  with  bars  are  the  corresponding 
original  matrices  and  vectors  after  the  usual  pivot  operations 
within  the  primal  and  dual  problems.  When  an  artificial  variable 
becomes  non-basic,  it  is  dropped  from  further  consideration. 
This  results  in  the  A  matrix  and  the  H  vector.   Each  iteration 
in  the  standard  primal-dual  algorithm  consists  of  the  following: 
a)   Finding  the  most  negative  d.  value  for  those  indices 
j  having  "c.  =0. 


b)  If  there  are  no  indices  with  this  property,  the 

multiplier  k  =  min  [c./  -  d.]  is  used  to  obtain  at 

d  .<0   J     3 
J 
least  one  c .  equal  to  zero  by 

c.  =  c.  +  k(d.) 
J    J      J 

c)  The  index  r  =  i  is  obtained  for  the  minimum  d.  with 

J 

c.  =  0.   Then  x   is  selected  to  be  basic. 

J  r 

d)  Then  min  (h./a.   )  is  found  for  a.    >  0.   This  will 

1   i,r  i,r 

select  a  pivot  row  in  either  the  X  or  Y  sets  of 

D 

equations . 

e)  The  simplex  pivot  method  is  then  used.   If  the  pivot 
row  is  in  the  Y  set  of  equations  the  artificial  variable 
in  that  row  is  dropped.   Thus,  in  any  case,  one  of  the 
forms  in  (2)  is  achieved. 

When  w  =0,  the  basic  solution  is  optimal  and  the  problem 
is  solved. 

The  actual  computations  should  include  rules  to  prevent 
cycling,  such  as  the  perturbation  technique  or  the  use  of 
lexicographic  ordering. 
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SECTION  III 

DESCRIPTION  OF  THE  NEW  TECHNIQUE 

The  primal-dual  algorithm  is  modified  in  a  two-step  procedure. 
In  the  first  step,  a  pivot  row  in  the  Y  set  of  equations  is  selected, 
By  only  considering  elimination  of  artificial  variables,  after  m 
or  less  iterations,  w  equals  zero.  At  this  point  the  equations 
have  the  form : 

X_  +  AX  =  H 

D 

z  +  CX  =  z 

o 

where  all  c.  are  positive.   If  all  the  h.  are  positive,  the  solution 
J  i 

is  optimal.  We  then  stop.   However,  in  most  cases,  this  stage  of 
the  problem  will  contain  one  or  more  negative  h..   In  this  case 
the  subscript  "it"  is  assigned  to  the  row  with  the  most  negative 
h..   The  following  procedure  is  then  used  as  the  second  step  of  the 
modification: 

a)  The  equation  containing  h   is  multiplied  through  by  minus 

one.   This  produces  h.   >  0. 

it 

b)  The  "it"  equation  is  then  added  to  all  equations  with 
negative  h. .  All  basic  variables  except  the  "it"  equation 
now  become  feasible. 

c)  To  maintain  the  canonical  form  an  artificial  variable  is 

added  to  the  "it"  equation.  A  new  infeasibility  equation 

with  w  =  h . ,_  is  formed  from  the  "it"  equation, 
o    it  n 

The  process  now  reverts  to  the  primal-dual  algorithm  to  achieve 
w  equal  to  zero.   Since  the  canonical  form  has  been  maintained  and 
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all  c  have  been  kept  positive  throughout  the  process,  once 
J 

w  equals  zero,  the  solution  is  optimal. 
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SECTION  IV 

COMPUTATIONAL  RESULTS 

Several  problems  were  run  on  an  IBM  360/67  using  both  the 
standard  primal-dual  and  the  primal-dual  with  the  new  technique, 
The  problems  used  were  test  problems  from  the  SHARE  Linear 
Programming  Project  and  two  highly  degenerate  assignment  type 
problems.   Results  are  given  in  Table  I  with  the  two  assignment 
problems  listed  first. 

TABLE  I 

TABLE  OF  COMPARATIVE  RESULTS 


PRIMAL- DUAL  WITH 


SIZE 


STANDARD 
PRIMAL- DUAL 


SIMPLEX  THE  NEW  TECHNIQUE 

PROBLEM  (m  x  n)        ITERATIONS      Iterations      Seconds      Iterations      Seconds 


Assignment  74  x  150 

Assignment  34  x  376 

SHARE    1A  33  x  64  7  0 

SHARE    ID  27  x  45  61 

SHARE    IE  31  x  106  66 

SHARE   IF  66  x  135  184 

SHARE   2A  30  x  103  114 


89 
74 
39 
38 
96 
125 
79 


32 

12 

7 

30 
86 
29 


83 

97 

39 

38 

136 

173 

119 


29 
28 

7 

42 

114 

43 


In   the   two  assignment   problems    the   primal-dual  with   the  new  technique 
resulted   in  optima  different   from   that   of  the   standard   primal-dual   since 


Philip  Wolfe   and  Leola  Cutler,    "Experiments    in  Linear   Programming," 
Recent  Advances   in  Mathematical   Programming    (New  York:     McGraw-Hill 
Book  Co.,    Inc.,    1963),    p.    198. 
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these  problems  have  non-unique  solutions.   The  two  methods  resulted 
in  the  same  optima  in  all  other  problems. 

No  accurate  timing  data  is  available.   The  timing  data  in  the 
above  table  though  does  give  an  indication  of  the  reduction  in 
execution  times  possible  with  the  new  technique.   However,  the 
modified  method  does  perform  less  work  than  the  usual  method 
because  in  the  first  step  ratios  and  comparisons  are  only  taken 
for  the  artificial  rows.   The  restart  of  the  problem  in  the  second 
step  requires  only  an  insignificant  amount  of  time.   Thus,  in  the 
table,  the  number  of  iterations  represents  a  measure  of  improvement 
of  the  new  method.   The  new  technique  shows  significant  reductions 
in  the  number  of  iterations  from  the  two-phase  standard  simplex 
method  in  four  out  of  the  five  problems  where  a  comparison  is 
possible.   When  comparing  the  new  technique  with  the  standard 
primal-dual,  a  similar  reduction  is  evident  in  the  larger  problems. 
It  therefore  appears  that  the  new  method  presents  a  significant 
improvement  when  the  number  of  variables  is  large. 
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SECTION  V 

DESCRIPTION  OF  THE  COMPUTER  PROGRAM 

The  accompanying  program  (see  Appendix)  is  written  in  the 

FORTRAN  IV  language,  and  was  used  to  run  all  test  problems  with 

the  exception  of  the  two  assignment  type  problems  (since  these 

only  deal  with  zeros  or  ones  in  the  A  matrix,  considerable 

reduction  in  the  core  memory  requirement  could  be  made  with  a 

slight  modification  of  the  program).   The  program  will  handle  up 

to  69  constraint  equations  and  up  to  149  variables,  and  requires 

about  160,000  bits  of  core  storage.   Columns  may  be  designated 

by  alphanumeric  names,  but  rows  must  be  consecutive  integer 

values  not  to  exceed  the  value  69.   Only  the  non-zero  elements 

of  A,  C,  and  H  need  be  read  in.  All  h.  must  be  positive,  but 

negative  c .  are  allowed. 
J 

Names  used  in  this  program  conform  to  those  already  outlined 
in  this  paper.   Those  not  previously  discussed  are  defined  as: 

FF  is  a  multiplicative  constant  used  to  insure  that  all  a.  . 

are  integer. 
EP  is  a  constant  used  in  the  round-off  error  control. 
DI  is  the  cumulative  multiplication  of  the  pivot  elements. 
CE  is  a  vector  identical  to  the  initial  C,  and  is  used  to 

calculate  z. 
AX  is  the  pivot  column  in  the  Revised  Simplex  tableau. 
CNAME  is  a  vector  of  column  names. 
DUAL  is  the  vector  of  dual  variables. 
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DB  is  a  vector  of  simplex  multipliers  used  to  calculate  the  d . . 
R  is  a  vector  used  to  record  whether  each  variable  is  basic 

or  artificial. 
E  is  a  computed  scalar  used  in  the  updating  and  round-off 

control. 
P  is  a  vector  of  multipliers  used  in  the  actual  pivot  operation. 
ITR  indicates  the  number  of  iterations  performed. 
INDEX  controls  the  use  of  the  new  technique  (INDEX  equals  zero 

yields    the   standard   primal-dual,    and   INDEX  equals   one 

incorporates  the  new  technique) . 
IB  is  a  vector  of  variable  designators  (a  negative  number 

indicates  an  artificial  variable,  and  a  positive  number 

indicates  a  basic  variable). 
All  other  names  are  defined  in  the  program. 

Once  the  data  has  been  read  in,  the  c.  are  tested  for  negative 

J 

values.   If  all  c.  are  positive,  the  primal-dual  commences.   However, 
if  there  is  one  or  more  negative  c.  the  following  procedure  is  fol- 


lowed: 


a)  RMIN  is  set  equal  to  the  minimum  c.. 

b)  n  is  set  equal  to  n  +  1  and  m  is  set  equal  to  m  +  1. 

c)  The  coefficients  of  all  variables  in  the  added  constraint 
equation  indicated  in  b)  above  are  set  equal  to  FF  including 
x  .   (Variable  x  does  not  appear  in  any  other  constraint 
equation.) 

d)  Then  c   is  set  equal  to  minus  RMIN. 

n 
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e)   The  right  hand  side  value,  h  ,  is  entered  as: 

m 

h  =  k(FF) 
m 

where  k  is  a  number  larger  than  the  sum  of  the  x. 

J 

in  the  optimal  tableau. 
The  problem  is  now  modified  so  it  may  be  handled  by  use  of  the 
primal-dual  algorithm. 

The  remainder  of  the  program  is  a  FORTRAN  representation  of 
the  steps  outlined  in  the  sections  "The  Standard  Primal-Dual 
Algorithm"  and  "Description  of  the  New  Technique."  However,  for 
computer  calculations,  an  essential  addition  to  the  procedures  is 
some  control  over  round-off  error. 

Round-off  error  develops  when  two  large  numbers,  say  y  and  z 
(differing  only  slightly  from  each  other),  are  subtracted  from 
each  other  such  as: 

x  =  y  -  z 
In  this  program  the  absolute  value  of  x  divided  by  y  is  compared 
to  EP.   If  it  is  less  than  EP,  x  is  set  equal  to  zero.   As  a 
further  control,  the  values  of  DB,  B,  and  H  are  updated  each 
iteration. 
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C      PROGRAM  TO  SOLVE  PRIMAL-DUAL  LINEAR 

C      PROGRAMMING  PROBLEMS 

C 

IMPLICIT    REAL*8    (A-H,0-Z) 

COMMON  A(70,15O),B(7O,70) ,C( 150) ,D( 150) ,H(70) ,AX( 7C) 

COMMON  0B(7C) , DUAL ( 70 ) , R ( 1 50 ) , CE ( 1 50 ) ,  E , W, Z , D I , FF , EP 

COMMON  CNAME( 150) 

COMMON  IB(70),M,N,JS,IR, INDEX , I TR , Ml , Nl 

FEAL*6  NINE/'9999     •/ 

FF=10000. 

EP=0.000C01 

DI  =  1. 

M=0 

DO  1  1=1,70 

H( I)=0. 

DO  1  J=1,15C 

1  A( I, J)=0. 

DO  2  1=1,150 
CNAME(I)=C. 
CE(I  )=C. 

2  C( I)=0. 
C 

C      READING  IN  COEFICIENTS  OF  CONSTRAINT  EQUATIONS 
C 

DO  5  N=l,1000 

READ  (5,3)  CB,  J,  RB 

IF(CB.EO.NINE)  GO  TO  16 

3  FORMAT  (6X,  A6,  16,  F10.4) 
00  4  1=1,150 

IF(CNAME(  I) .EQ.CB)  GO  TO  30 
IF(CNAME(  I) .EO.O. )  GO  TO  30 
GO  TO  4 
30  A(J,I)=RB*FF 
CNAMECI )=CB 

WRITE  (6,102)  I,  J,  A(J,I) 
102  FORMAT  (IX,  2110,  F18.9) 
IF(J.LE.M)  GO  TO  5 
M=J 
GO  TO  5 

4  CONTINUE 

5  CONTINUE 
C 

C      READING  IN  RIGHT  HAND  SIDE  VALUES 
C 

16  DO  6  N=1,1C0 
READ  (5,14)  I,  RB 

14  FORMAT  (  118,  F10.4) 
IFU.GT.999)  GO  TO  17 
H(I)=RB*FF 
WRITE  (6,14)  I,  H(I) 

6  CONTINUE 
C 

C      REAOING  IN  COSTS 
C 

17  DO  9  J=1,20C 
PEAD  (5,7)  CC,  II,  RB 

">,  F10.4) 
TO  18 


f tAU  \ D , f I   Lt»   11, 

7  FORMAT  (6X,  A6,  16, 
IF(CC.EQ.NINE)  GO  T 
DO  11  K=l,150 
IF(CC.EQ.CNAME(K) ) 
GO  TO  11 


GO  TO  12 
GO  TO  11 
12  C(K)=RB*FF 
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II,  C(K) 
F15.5) 


CE(K)=RB 

WRITE  (6,13)  K, 
13  FORMAT  (6X,  216, 

GO  TO  9 
11  CONTINUE 
9  CONTINUE 
18  CONTINUE 

DO  20  1=1,150 

IF<CNAME<  U.GT.O.)    GO    TO    20 

N=I-1 

GO  TO  21 

20  CONTINUE 

21  WRITE  (6,101)  M,  N 

101  FORMAT  (10X,  'NUMBER  OF  ROWS=',  15,  15X, 
l'NUMBER  OF  COLUMNS=»,  15) 
WRITE  (6,103)  (CNAME(I),  1=1,175) 
FORMAT  (IX,  20A6) 


103 
104 


rurvriMi       »  i.  a  •       cuhui 

WRITE  (6,104)  (C(J),  J=1,N) 
FORMAT  (IX,  10E12.4) 


INDEX=1 

N1=N+1 

M1=M+1 

CALL  LIP 

STOP 

END 
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c 

c 


SUBROUTINE  LIP 


IMPLICIT  RFAL*8  <A-HfO-Z) 

COMMON  A<70,150),B(70,70),C(150),D(150) ,H<7C) ,AX(70) 

COMMON  DB(7C) , DUAL ( 70 ) , R( 150 ) ,CE<150),F,W,Z,DI,FF,EP 

COMMON  CNAME( 150) 

COMMON  IB<70),M,N,JS, IR , INDEX, I TR, M1,N1 

E=0.5/DI 
C 

C      MODIFYING  THE  PROBLEM  IF  THE  OBJECTIVE 
C      FUNCTION  CONTAINS  NEGATIVE  COSTS 
C 

RMIN=0. 

DO  50  1=1, N 

IF(C(  D.GE.E)  GO  TO  50 

IF(C<I).GE.RMIN)  GO  TO  50 

RMIN=C( I) 
50  CONTINUE 

IF(RMIN.EO.O. )  GO  TO  60 

DO  52  1=1, N 

52  C(I  )=C(  D-RMIN 
N=N+1 

M  =  M+1 

H(M)=1000.*FF 
C(N)=-RMIN 
DO  53  J  =  1,N 

53  A(M,J)=FF 
C 

C      INITIALIZING  THE  TABLEAU 
C 

60  DO  3  J=1,M 
DO  2  1=1, M 

2  B(I,J)=0. 

3  B(J,J)=1. 
DO  4  1=1, M 
DB(I )=-l. 
DUAL( I»=0. 

4  £B(I)=-I 
DO  5  J=1,N 
R( J)=C. 

0( J)=C. 

DO  5  1  =  1,  M 

D( J)=D(J)-A( I, J) 

5  CONTINUE 
W=0. 

DO  6  1=1, M 

6  W=W+H( I) 
ITR=0 

C 

C      OBTAINING  A  COST  EQUAL  TO  ZERO  AND  A  PIVOT  COLUMN 

C 

7  IF(DABS(W)-E)100,100,8 
3  DO  9  J=l,N 

IF(D( J) )10,9,9 
9  CONTINUE 
GO  TO  101 
10  JS=0 

00  11  J=1,N 

IF(R(J).E0.1.)  GO  TO  11 
IF(DU).GE.-E)  GO  TO  11 
IFUS.EQ.O)  GO  TO  13 
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12  IF<CX*C(J)/D(  J)  )11,U,13 

13  CX=-C(J)/D( J) 
JS=J 

11  CONTINUE 
47  DO  14  J=1,N 

14  C< J)=C( J)+CX*0(J) 
C( JS)=0. 

00  301  J=1,N 
IF(C(J).LT.E»  GO  TO  301 
GO  TO  302 

301  CONTINUE 
GO  TO  100 

302  Z=0. 

DO  303  1=1, M 

J=IB(I) 

IF(J.LE.C)  GO  TO  303 

Z=Z+CE( J)*H( I  J 

303  CONTINUE 

DO  40  1=1, M 

40  DUAL( I)=DUAL< I )-CX*DB( I) 

15  IF(W-E)100,100,41 

41  DO  16  J=1,N 

1  F(C( J) -0.00001) 17, 17, 16 

17  C(J)=0. 

18  IF(D( J)-O(JS) )19,16,16 

19  JS=J 

16  CONTINUE 

WRITE  (6,213)  JS,  D(JS) 
213  FORMAT  (1H,  »JS=',  14,  3X,  • D= • ,  F12.4) 
IF(D(JS) )2C,7,7 

20  CALL  ELEM 
IF(DI-1.I528, 26,26 

26  CALL  PIVOT 
GO  TO  15 
C 

C      PRINTING  OF  RESULTS 
C 

100  WRITE  (6,204) 

204  FORMATdH  ,  'RESULT') 
GO  TO  425 

101  WRITE  (6,205) 

205  FORMATdH  ,  •  INFE  AS  I8LE  •  ) 
STOP 

425  DO  426  1=1, M 

426  WRITE  (6,526)  IB(I),  H(I),  DUAL(I) 

526  FORMAT  (1H  ,  'BASIS  VARIABLE  NUMBER=»,  14,  5X, 
l'VALUE  OF  BASIS  VARIABLE=',  F12.4,  5X, 
l'DUAL=',  F12.4) 

IF(RMIN.EO.C)  GO  TO  427 
ZZ=0. 

DO  429  K=1,M 
J=IB(K) 

IF(J.LE.O)  GO  TO  429 
ZZ=ZZ+CE< J)*H(K) 
429  CONTINUE 
GO  TO  428 

427  ZZ=Z 

428  WRITE  (6,527)  W,  ZZ 

527  FORMAT  (1H  ,  • W= • ,  F12.8,  5X,  »Z=',  F12.5) 
IFdNDEX.EO.il  GO  TO  1 

528  RETURN 
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c 

C      THE  NEW  TECHNIQUE 


C 


1  IT  =  0 

DO  30  1=1, M 

IF(H(  D.GT.-E)  GO  TO  30 

IF(IT.EO.O)  GO  TO  31 

IF(H( D.GE.H< IT) )  GO  TO  30 
31  IT  =  I 
30  CONTINUE 

IF(IT.EO.O)  GO  TO  528 

IttC IT>=-IT 

H(IT)=-H( IT) 

DO  34  K=1,M 

DB(K)=B<  IT,K) 

34  B( IT,K)=-B< IT,K) 
DO  33  1=1, M 

IF(H(  D.GT.-E)  GO  TO  33 
H(D=H(  D+H(  IT) 
DO  35  K=1,M 

35  B(I,K)=B<  I,K)+B<  IT,K) 
33  CONTINUE 

W=H( IT) 
DO  36  J=i,N 
AX(IT)=0. 
DO  21  K=1,M 
21  AX( IT)  =  AXUT)+B( IT,K)*A(K, J) 

36  D( J)=-AX( IT) 
INDEX=0 

GO  TO  8 
END 
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c 

c 


SUBROUTINE  ELEM 


IMPLICIT  REAL*8  (A-HfO-Z) 

COMMON  A(  70,  1  50  ».  R  (  70,  70)  ,C(  150.0(150)  ,H(70)  ,AX(70) 

C OMMON  DB  (  7  0 )  ,  DU  AL  (  70  )  ,  R  (  1  50  )  ,  C  E  (  1  50  )  ,  E  ,  W  f  Z  ,  D I  ,  F  F  ,  E  P 

COMMON  CNAME( 150) 

COMMON  IB (70) ,M,N,JS, I R , INOEX , I TR , Ml , Nl 

IR=0 

J  =  JS 
C 

C      CALCULATING  THE  REVISED  SIMPLEX  TABLEAU  PIVOT  COLUMN 
C      WITH  ROUNO-OFF  ERROR  CONTROL 
C 


DO  2  2  I=1,M 

0M=0. 

DP  =  0. 

AX(I )=0. 

DO  21  K=1,M 

T=B( I,K)*A(K,J) 

IF(T.GE.O. )  GO  TO  40 

DM=DM+T 

GO  TO  21 

40 

DP=DP«-T 

21 

CONTINUE 

AX(I )=DP+DM 

IFCOP.EO.O. )  GO  TO  22 

IF(OABS(AX( I)/DP).GT.EP> 

AX(I )=0. 

22 

CONTINUE 

GO  TO  2  2 


C 

C      FINDING  THE  PIVOT  ROW 
C 


GO  TO  25 
70  IF( IB(I))74,25,73 

73  IF(IB(IR))25,25,4 

74  IF(IBUR)  )4,25,3 
4  IF(AX( IR)-AX(  I  ))25,25, 

25  CONTINUE 


25  CONTINUE 
RETURN 
END 
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SUBROUTINE  PIVOT 
C 
C 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMON  A  (70,150,  8  (  70,70  ),C<  150,0(150)  ,H(  70)  ,AX(  70) 

COMMON  DB(70) , DUAL ( 70 ) , R ( 1 50 ) , C E \ 1 50 ) , E , W, Z , DI , FF  ,  EP 

COMMON  CNAMEC150) 

COMMON  IB(70),M.N,JS,IR, INOEX, ITR, M1,N1 

OIMENSION  P(200) ,BX(200) 

WRITE  (6,214)  IR,  AX(IR) 
214  FORMAT  ( 1H,  25X,  MR=«,  13,  3X,  »AX=',  F12.4) 
C 

C      PERFORMANCE  OF  THE  PIVOT 
C 

43  IFUB(IR))  1,1,2 

2  K=IB(IR) 
R(K)=C. 
GO  TO  3 

1  P(JS)=1. 

3  IB(IR)=JS 
DI=DI*DABS( AX< IR)  ) 
IFtOI.GE. 9999998.)  GO  TO  40 
DI  =  IDINT(DI-»-0.5) 

40  CONTINUE 
E=0.5/DI 
DO  27  1=1, M 

27  F(  I)=-AX( I)/AX(IR) 
P( IR)=1./AX(IR)-1. 
P(M+1)=-D(JS)/AX( IR) 
DO  28  J  =  1,M 
BX(J)=B( IR,J) 

28  CONTINUE 
HX  =  H(  IR) 

DO  29  J=1,M 

DO  29  1  =  1, M 

T=B( I, J) 

Bdt  J)=B(  I,J)+P(I)*BX(J) 

IF(T.EO.O.)  GO  TO  29 

IF(DABS(B(I,J)/T).GT.EP)  GO  TO  29 

B( I, J)=0. 

29  CONTINUE 

DO  30  1  =  1, M 

T^H(  I) 

H(I)=H( I)+P(I)*HX 

IF(T.EO.O.)  GO  TO  30 

IF(DABS(H(  I)/T).GT.EP)  GO  TO  30 

H(I)=0. 

30  CONTINUE 
W=W-P(M*1)*HX 
DO  31  J=1,M 
T=DB( J) 

DB( J)=DB(J)^P(M+1)*BX(J) 
IF(T.EO.O.)  GO  TO  31 
IF(OABS(DB( J)/T).GT.EP)  GO  TO  31 
DB( J)=0. 

31  CONTINUE 
C 

C      UPDATING  THE  DB,  B,  AND  H  VALUES 


C 


IF(DI.GE. 999998.)  GO  TO  81 

DO  92  J=1,M 

DF  =  DI*DABS(DB(JM+0.5 
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IF(OF.GT. 999998.)  GO  TO  61 
OF=IDINT(DF)/DI 
D8( J)=OSIGN(DF,DB(J)  I 
61  00  92  1=1, M 

DF=DI*DABS(B<  I, J)  )+0.5 

IF(DF. GT. 999998. )  GO  TO  92 

DF=IDINT(DF)/DI 

B(I, J)=DSIGN(OF,B<  I, J)  ) 

92  CONTINUE 

00  93  1=1, M 

DF=DI*DABS(H(  I  )  )  +0 . 5 

IF(DF.GT. 999998.)  GO  TO  93 

CF=IDINT(OF)/OI 

H( I)=DSIGN(OF,H( I )  ) 

93  CONTINUE 
81  CONTINUE 

DF=W*DI 

IF(OF. GE. 9999998. )  GO  TO  87 

W=IDINT(DI*W+0.5)/DI 
87  CONTINUE 
C 
C      CALCULATION  OF  D(J)»S  FOR  NEXT  ITERATION 


C 


DO  34  J=1,N 
DM=0. 
0P=0. 
D( J)=0. 
DO  33  K=1,M 
T=DB(K)*A(K,J) 
IF(T.GE.O.)  GO  TO  41 
OM=DM+T 
GO  TO  33 
41  DP=DP+T 

33  CONTINUE 
D< J)=DP+DM 

IF(DP.EO.O. )  GO  TO  34 
IF(DABS(D(J)/DP).LE.EP)  GO  TO  35 
IF(DABS(D(J) )-E)35,35,34 

35  D(J)=0. 

34  CONTINUE 
D(JS)=0. 

IFCW.GT. 0.001)  GO  TO  44 

w=o 

44  ITR=ITR+1 

WRITE  (6,210)  ITR,  W,  Z,  01 
210  FORMAT  ( 1H  ,  5X,  MTR=«,  15,  5X,  • W= • ,  F18.8 
15X,  «Z=»,  F12.4,  5X,  '01=',  E12.4) 

20  RETURN 

21  STOP 
FND 
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